Functional imaging through scattering medium via fluorescence speckle demixing and localization

Recently, fluorescence-based optical techniques have emerged as a powerful tool to probe information in the mammalian brain. However, tissue heterogeneities prevent clear imaging of deep neuron bodies due to light scattering. While several up-to-date approaches based on ballistic light allow to retrieve information at shallow depths inside the brain, non-invasive localization and functional imaging at depth still remains a challenge. It was recently shown that functional signals from time-varying fluorescent emitters located behind scattering samples could be retrieved by using a matrix factorization algorithm. Here we show that the seemingly information-less, low-contrast fluorescent speckle patterns recovered by the algorithm can be used to locate each individual emitter, even in the presence of background fluorescence. We test our approach by imaging the temporal activity of large groups of fluorescent sources behind different scattering phantoms mimicking biological tissues, and through a brain slice with a thickness of ~ 200 μm.

: Rank estimation from experimental data. For the same recorded dataset, we show the average residual error for five different NMF realizations (with different random initializations) for different rank values. After the rank surpasses the true number of sources in the sample, the residual error decreases at a much lower rate, a phenomena that can be used to estimate the number of emitters in the sample.

NMF inversion problem
In order to solve the general NMF problem, multiple numerical methods can be used [34,35]. While many of the currently available solvers simply tackle the simplest form of the inversion problem found in the main text (min I − W · H 2 F subject to W, H > 0), it is possible to add regularization terms with some a priori information about the system, such as the sparsity of either the fingerprints or the temporal activities of the sources in the sample. Thus, we can formulate the NMF problem as: min W,H>0 0.5 · ||I − W H|| 2 β + α W · l 1r · n pixels · ||vec(W )|| 1 + α H · l 1r · n f rames · ||vec(H)|| 1 where ||A|| 2 F = i,j A 2 ij corresponds to the Frobenius norm of a matrix, A, ||vec(A)|| 1 = i,j abs(A ij ) corresponds to the l 1 element-wise norm, and ||I − W H|| β represents the desired β-norm to calculate (1, 2, ≤ 0). Here, the additional terms introduced account for both the sparsity of W and H, and are governed by α W , α H , and l 1r (with two scaling factors, n pixels and n f rames , accounting for the vast differences between the number of elements of W and H). Both α W and α H can take different values in order to weight the strength of the regularization between W and H, and l 1r can be used to continuously choose between different penalty forms. For the limit l 1r = 0, the penalty behaves like a standard Frobenius norm, while l 1r = 1 corresponds to an element-wise l 1 penalty (favoring sparsity). In our case, we acquire low contrast images that result of the incoherent addition of many highly-contrasted individual speckles. Furthermore, these patterns do not fully cover the field-of-view of the camera, so some degree of sparsity is to be expected on each individual fingerprint. Moreover, this promotes recoveries where the fingerprints have higher contrast, which greatly helps the localization procedure. Last, the temporal activities that we use to mimic neuronal activity consist of short bursts of activity, usually followed by longer decay times and periods of little activity, so it is reasonable to consider some sparsity on the recovery of H. We find that, in our experimental conditions, a good compromise between fidelity, sparsity, and recovery time is found with β = 2, α W = 1.5, α H = 0.5, and l 1r = 0.5. While these regularization parameters have to be manually tuned and are experiment-dependent, several approaches to automatically estimate their values could be explored in the future [36,37]. The full system is solved by using the scikit-learn NMF package [29], and it takes a few minutes to compute for datasets consisting of 500 frames with resolutions in the order of 300 × 300 pixels using a desktop CPU (Intel i7-9700) with 16 Gb of RAM. Bigger datasets and/or faster reconstruction times could be reached by using GPU-based implementations, but this lies outside of the scope of this work.
Step-by-step localization procedure Here, we introduce the post-processing workflow to obtain the location of the emitters from the recorded dataset. First, we crop and filter the frames recorded by the camera (as the sensor is larger than the area covered by the speckle patterns). Then, we remove the intensity envelope by high-pass filtering. In the experiments where there is a constant signal present (as in Fig.3), we perform a rank-1 NMF to identify this constant component in the dataset, which we later use to initialize both W and H when performing a full-rank NMF with the rank set to the estimated number of emitters. This helps unmix the time-varying fingerprints and the background present in all the frames due to the constant fluorescence signal. Otherwise, we just initialise the NMF with the Nonnegative Singular Value Decomposition (NNSVD) of the recorded dataset. After the NMF is performed, we deconvolve all the fingerprints in pairs to locate the shifts between them, and finally we merge all the information in the full location map. The codes can be found at [31].

Algorithm 1:
Step-by-step localization procedure Result: Returns the location map of the emitters in the sample by performing NMF over the recorded dataset.
The fingerprints provided by the NMF are deconvolved to calculate the distances between the different sources.
Post-process recorded dataset: select the region of the sensor with speckle patterns (cropping), perform binning (reduce size to increase speed) and high-pass filtering (remove envelope, increase contrast) if constant background = True then Do Rank-1 NMF to estimate constant background in the dataset end Perform a full-rank NMF on the recorded dataset (n s = number of sources) if constant background = True then Set rank = n s + 1 Perform the NMF, initializing W and H with the result of the rank-1 NMF else Set rank = n s Perform the NMF, initializing W and H with the Nonnegative Singular Value Decomposition (NNSVD) of the recorded dataset end Calculate source positions by deconvolving the fingerprints provided by the NMF algorithm. for i = 1 : n s do for j = 1 : n s do Deconvolve f ingerprint i and f ingerprint j δ(x − x i,j 0 , y − y i,j 0 ) = w i * −1 w j Calculate distance between source i and source j as the position of the delta-like peak (x i,j 0 , y i,j 0 ) end Combine deconvolutions in a partial location map (M i ): Correct shifts between all the partial location maps by using the distances between the sources, then add together to generate the full location map